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ABSTRACT. The "equation-free" approach has been proposed in recent years as a general 
framework for developing multiscale methods to efficiently capture the macroscale behav- 
ior of a system using only the microscale models. In this paper, we take a close look at 
some of the algorithms proposed under the "equation-free" umbrella, the projective inte- 
grators and the patch dynamics. We discuss some very simple examples in the context 
of the "equation-free" approach. These examples seem to indicate that while its general 
philosophy is quite attractive and indeed similar to many other approaches in concurrent 
multiscale modeling, there are severe limitations to the specific implementation proposed 
by the equation-free approach. 



1. Introduction 

The purpose of this note is to examine some of the basic issues surrounding the "equa- 
tion-free" approach proposed in [17], which has been pursued in recent years as a general 
tool for multiscale, multi-physics modeling. To begin with, the equation-free approach 
is an example of concurrent coupling technique. In contrast to sequential coupling tech- 
niques which require establishing the macroscale equations through precomputing, con- 
current coupling techniques compute the required macroscale quantities "on-the-fly" from 
microscopic models [1, 2]. The most well-known example of such concurrent coupling 
techniques is perhaps the Car-Parrinello molecular dynamics which computes the atomic 
interaction forces "on-the-fly" by solving the electronic structure problem [5]. Other al- 
gorithms, such as the extended multi-grid method [4] and the heterogeneous multiscale 
method (HMM) [8] are all example of the concurrent coupling approach. 

At a technical level, a key idea in the "equation-free" approach is to make use of scale 
separation in the system. There are many different ways of exploiting scale separation. 
In [6] and [5], time scale separation was used to artificially slow down the time scale of 
the microscopic system. As for spatial scales, homogenization-based methods (such as 
the ones that use representative averaging volumes [3]) and the quasicontinuum methods 
[16] are all examples of algorithms that explore the separation of spatial scales. Most 
closely related to the "equation-free" approach is perhaps the extended multi-grid method 
[4]. In his review article [4], Ac hi Brandt described ideas that can be used to extend 
multi-grid techniques to deal with multiscale, multi-physics problems in order to capture 
the macroscale behavior of a system using microscopic models such as molecular dynam- 
ics. As is common in multi-grid methods, the ideas of Brandt rely heavily on mapping 
back and forth between the macro- and micro-states of the system, through prolongation 
and restriction operators (which are called respectively reconstruction and compression 
operators in HMM, and lifting and restriction operators in the "equation-free" approach). 
Brandt realized that central to the efficiency of these algorithms is the possibility of only 
performing microscopic simulations in small samples for short periods of times, as a result 
of the scale separation in the system. These ideas were later adopted by both HMM and 
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the "equation-free" approach. In fact, HMM and the "equation-free" approach are both 
alternative approaches with the same motivation and objective. 
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micro to Macro 


Extended multi-grid 


interpolation 


restriction 


HMM 


reconstruction 


compression 
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restriction 



While the general philosophy of the "equation-free" approach is very similar to the ex- 
tended multi-grid method and HMM, the "equation-free" approach proposes its own ways 
of implementing such a philosophy, in particular, ways of dealing with scale separation. 
The basic idea is to use extrapolation in time and interpolation in space. More precisely, 
two important building blocks of the "equation-free" approach are: 

i) Projective integrators: (An ensemble of) the microscale problems are solved for a 
short period of time using small time steps. The time derivative of the macro vari- 
able is computed from the results of the last few steps and then used to advance 
the macro variable over a macro time step. It is easy to see that such a procedure 
amounts to extrapolation, and indeed the authors state in [13]: "The reader might 
think that these should be called 'extrapolation methods,' but that name has al- 
ready been used [...]. Hence we call the proposed methods projective integration 
methods." 

ii) The gap-tooth scheme: The microscopic problem is solved in small domains (the 
teeth) separated by large gaps. The solution is averaged over each domain and 
then interpolated to give the prediction over the gaps. 

The combination of these two ideas gives directly the so-called "patch dynamics" [17]. 

Detailed understanding of the "equation-free" algorithms is made difficult by the fact 
that the "equation-free" papers are generally quite vague. The present note should be re- 
garded as an attempt to pin down some of these details. Indeed this was initially intended 
as a regular journal article. But it soon becomes clear that there is still substantial disagree- 
ment between our understanding of the "equation-free" approach and that of its developers. 
However, we believe the simple examples that we discuss here do shed some light on the 
"equation-free" approach and should be made available to a larger audience in some form. 
We are grateful to Yannis Kevrekidis for a detailed report on the earlier version of this 
note. Some of his comments have been taken into account in this revised version. We 
also welcome any discussion about the issues raised in this note, the most important of 
which being: What really is the "equation-free" approach? Indeed our primary purpose of 
presenting this note is to prompt such a discussion. 

2. Projective integrators for stochastic ODEs 

Projective integrators were proposed as a way of extrapolating the solution of an explicit 
ODE solver for systems with multiple time scales using large time steps. The basic idea 
is to run the microscopic solver (using small time steps) for a number of steps, and then 
estimate the time derivative and use that to extrapolate the solution over a large time step 
[13]. For stiff ODEs, the extrapolation step is applied to the whole system [13]. For general 
multiscale problems, the extrapolation step is applied only to the slow variables [17, 15]. 

In the case of stiff ODEs, projective integrators can give rise to useful numerical sche- 
mes, as was demonstrated in [13]. In this case, the idea becomes very close to the ones 
proposed by Eriksson et. al for developing explicit stiff ODE solvers [12]. The objectives 
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of the two papers are quite different: For Eriksson et al, the objective is to find explicit and 
efficient stiff ODE solvers. For Gear et ai, the objective is to deal with general multiscale, 
multi-physics problems. However, in the general case such as the case considered in [15], 
projective integrators have serious limitations, as we now show. 

Denote by x the coarse variable of the system. The coarse projective integrator proposed 
in [15] performs the following steps at each macro time step (of size At): 

i) Create an ensemble of N microscopic initial conditions consistent with the known 
coarse variable x n at time step n. 

ii) Run the microscopic solver with these initial conditions for a number of steps, 
say k, with time step St. Denote the corresponding values of the coarse variables 
as Xj (kSt) where j = 1, . . . , N. 

iii) Perform ensemble averaging to get an approximation to the coarse variable. For 
example, 

1 N 

(1) x = -j2*^ k6 v 

3=1 

iv) Use this value to extrapolate the coarse variable to a time step of size At. 

(2) x n+1 =x n + At X ~ ' 



kSt 

Now consider the simple case when the coarse variable obeys effectively a stochastic 
ODE: 

(3) dx(t) = b(x(t))dt + dW{t) 
Since, to O(kSt), we have 

(4) xj (kSt) - x n = kSt b(x n ) + VkSt 

where {£"}, j = 1, ■ ■ ■ ,N are N independent Gaussian variables with mean and vari- 
ance 1, (2) becomes, to leading order 

At 1 N 

(5) x n+l =x n + Atb{x n )+ 



Mi N-. , 

(5) is equivalent in law to 

(6) x n+1 =x n + Atb(x n ) + —^=C, 

y/7jkSt 

where £™ is a Gaussian variable with mean and variance 1 . 

It is obvious from this that the effective dynamics produced by the coarse projective 
integrator depends on the numerical parameters N, k, 5t, and At. In particular, if NkSt ^> 
At, then the noise term in (3) is lost in the limit. If NkSt <C At, then the noise term 
overwhelms the drift term. In either case, one obtains a wrong prediction for the effective 
dynamics of the coarse variable. 

The only way to get a scheme consistent with (3) is to choose the numerical parameters 
so that they precisely satisfy NkSt = At. The reader should be aware, however, that this 
choice is not advocated in [15] and is, in fact, quite orthogonal to the original equation-free 
philosophy since it requires knowing beforehand that (3) is an SODE and not something 
else. In addition, it is easy to see that using NkSt = At leads to no gain in efficiency: The 
total cost is comparable to solving the microscopic problem in a brute force fashion using 
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St as the time step, since the size of the ensemble is equal to the number of microscopic 
simulation time intervals during a time duration of At: N = At/(kSt). 

For the case when NkSt » At, one might think of using the coarse projective integra- 
tors (or coarse molecular dynamics) as a way of simulating the dynamics dx/dt = b(x) 
in the context of molecular dynamics simulations. In this case the unknown drift b(x) is 
related to the gradient of the free energy and simulating dx/dt = b(x) is then a way to 
explore this free energy. Indeed, this appears to be how the scheme was actually used in 
[15]. The problem, however, is that using NkSt At leads again to a scheme which is 
no less expensive than a brute force solution of N replica of (3) using St as the time step. 

The problem above seems to be intrinsic to projective integrators in the context of SDEs 
because it is inherent to the fact that the dynamics (3) is dominated by the noise on short 
time scales and the extrapolation step in the projective integrators amplifies these fluctua- 
tions. Averaging them out can only be done at a cost which is at least comparable to the 
cost of a direct scheme. 



3. Patch dynamics 

Patch dynamics is proposed as a way of analyzing the macroscopic dynamics of a 
system using microscopic models. Like the extended multi-grid methods [4] and HMM 
[8, 11], it is formulated in such a way that scale separation can be exploited to reduce 
computational cost. 

The setup is as follows. We have a macroscale grid over the computational domain. The 
grid size Ax is chosen to resolve the macroscale variations but not the microscale features 
in the problem. Each grid point is surrounded by a small domain (the "tooth"), the size of 
which (denoted by h) should be large enough to sample the local microscale variations but 
can be much smaller than the macroscale grid size if the macro and microscales are very 
much separated. 

Given a set of macroscale values at the macroscale grid points, {U™}, at the n-th time 
step t n — nAt where At is the size of the macroscale time step, patch dynamics computes 
the update of these values at the next macroscale time step, {U™ +1 }, using the following 
procedure: 

i) Lifting: From {U™}, reconstruct a consistent microscopic initial data, denoted by 
u . 

ii) Evolution: Solve the original microscopic model with this initial data uo over the 
small domains (the "teeth") for some time St: ust = SstUQ. 

iii) Restriction: Average the microscale solution u$ t over the small domains. The 
results are denoted by {Ug t }- 

iv) Extrapolation: Compute the approximate derivative and use it to predict {U™ +1 }: 

fjn _ Tjn 

(7) U n+1 = U n + At--^ — — 



St 



or more generally: 



77™ — 77" 

(8) U n+1 = U n + At-^ — 

(1 — a)8t 

where a is some numerical parameter between and 1 . 

There are very few examples of how to implement these steps in practice. [17, 25] 
suggest the following: 
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For the lifting operator, in the small domain around the macro grid point Xj, use the 
approximate Taylor expansion: 

d 1 

(9) u (x) =^-£ fc (x-x,) fe 

k=0 

Here Dk is some approximations to the derivatives of the macroscale profile at Xj, for 
example: 

(10) ° 2 ' Dl 2A^ ' D ° U > 24 k D2 

Below we will consider the case when d = 2. 

When solving the microscale problem, [25] suggest extending the microscale domain to 
include some buffer regions in the hope that this would allow the use of any boundary con- 
ditions for the microscopic solver: By choosing sufficiently large buffers, the effect would 
be as if the microscale problem is solved in the whole space where and when averaging is 
performed. This introduces another spatial scale H which is the real size of the region on 
which microscale problems are solved. (The parameter h, which is (much) smaller than H, 
is the size of the domain over which averaging is performed). In the following discussion, 
we will take H to be infinity. 

Let us now examine this algorithm in more detail, using some very simple examples. 
Let us first consider the heat equation 

(11) d t u = d 2 x u 

For simplicity, let us assume Xj — 0. Denote u = D + D\X + \D 2 x 2 . We have 

(12) SstMx) = D + D lX + D 2 (^x 2 + St) 

Denote by Ah the averaging operator over the small domain (of size h), we have 

(13) U2t - A h S St u (x) = D + D 2 5t + ^D 2 h 2 = U n + D 2 5t 
Inserting this expression in (7) gives the familiar scheme: 

(14) U n+1 = U n + AtD 2 

as was shown in [17]. This is both stable and consistent with the heat equation, which is 
the right effective model at the large scale. 
Now let us turn to the advection equation 

(15) d t u + d x u = 
In this case, we have 

(16) S st u {x) =00 + 0^- St) + \v 2 {x - St) 2 
Hence, 

(17) = A h S St uo{x) =Do- D 1 St + X -D 2 St 2 + ^D 2 h 2 =U n - D x St + X -B 2 St 2 
and (7) becomes 

(18) U n+1 = U n + At(-Di + ^StD 2 ) 
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Since St <C At, the last term is much smaller than the other terms, and we are left essen- 
tially with a scheme which is unstable under the standard CFL condition that At <~ Ax: 

(19) U n+1 = U n - AtD 1 

due to the central character of D\. 

Aside from the stability issue, there can also be problems with consistency. Consider 
the following example: 

(20) d t u= -d\u 

The macroscale model is obviously the same model. However, it is easy to see that if we 
follow the patch dynamics procedure with d = 2, we would be solving dtU = 0, which is 
obviously inconsistent with the correct macroscale model. 

For these simple examples, the difficulties discussed above can be fixed by using dif- 
ferent reinitialization procedure for the micro-solvers. For the example of the convection 
equation, one should use one-sided interpolation schemes in the spirit of upwind schemes. 
For the last example, one should use piecewise 4th order polynomial interpolation. But 
in general, finding such a reinitialization procedure seems to be quite a daunting task, 
since it depends on the nature of the unknown effective macroscale model. Imagine that 
the microscopic solver is molecular dynamics. The reinitialization procedure has to take 
into account not only consistency with the local macrostates of the system (which is the 
only requirement for the extended multi-grid method and HMM), but also the effective 
macroscale scale model (which is unknown) such as: 

i) The order of the macroscale equation. 

ii) The direction of the wind, if the effective macroscale equation turns out to be a 
first order PDE. 

iii) Other unforeseeable factors. 

Indeed it is not at all clear how patch dynamics would work if molecular dynamics models 
are used to model macroscopic gas dynamics. 

To overcome these problems, the "equation-free" developers propose to design a num- 
ber of numerical tests to find out the nature of the effective macroscale equations. One such 
a procedure, the "baby-bathwater scheme" will be discussed in the next section. However, 
in addition to the technical issues, it is not clear what set of characteristics that we are 
supposed to test on. 

4. The "baby-bathwater scheme" 

As the last example shows, it is useful at times to know the order of the highest order 
derivatives that appear in the effective macroscale equation, even if we do not know all the 
details of the macroscale model. An algorithm was proposed in [19] for this purpose. The 
"baby-bathwater scheme", as it was called, promises to find the highest order derivative in 
the effective macroscale model, by performing simulations using the microscopic model: 
Assume that the effective macroscale model is of the form 

(21) d t U = F{U,d x U,--- ,d™U) 

The objective is to find m. 

This problem can be formulated abstractly as follows. Assume we have a function 
F = F(xi,X2, • • • ) and we know that it only depends on finitely many variables: F = 
F(x\,X2 ••• , x m ). Assume that we can evaluate F at any given point, can we find the 
value of m efficiently? 
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The basis idea used in [19] is that if F depends truly on the variable Xj, then the variance 
of F as Xj changes should not vanish. The practical difficulty is how to use this idea 
efficiently. 

Without getting into the details of the algorithm presented in [19], it seems quite clear 
that there is no fool-proof inductive procedure for finding m. Take an extreme case, say, 
F = F(xi, xioo). Without having the prior knowledge that F may depend on Xioo, an 
inductive procedure would likely conclude that F is only a function of x\. 

This example is of course quite extreme, and most practical situations are not like this. 
Nevertheless, it does raise some questions about the robustness of the algorithm presented 
in [19]. There is a much more serious concern, and this is associated with the well-known 
phenomenon that the order of the effective macroscale model depends on the scale we 
look at. A physically intuitive example is convection and diffusion of tracer particles in 
the Rayleigh-Bernard cells: At the scale of the cells, the tracer particles are convected and 
diffusion can be neglected. Hence the effective model is a first order equation. At scales 
much larger than the size of the cells, diffusion dominates. Hence the effective model is a 
second order equation. This means that the output of the "baby-bathwater scheme" should 
depend on the numerical parameter A in the scheme. 

This behavior can be demonstrated rigorously using the well-known results of Kesten 
and Papanicolaou [20]. Consider the dynamics if inertial particles in a stationary random 
force field: 

d 2 x 

(22) . FW 
In phase space, we can write this as 

( dx 

(23) ^ f 

dv cv ^ 

The density of the particles (in phase space) obeys the Liouville equation: 

(24) d t p + vd xP + F(x)d vP = 
However, if we consider the rescaled fields: 

{dx 5 _ I s 



(25) 



it was shown in [20] under some conditions on the random field F that as 5 — > 0, the 
process v s (-) converges to a diffusion process. In other words, if we consider the density 
of the particles in v space, then in this scaling we have 

(26) d t p + d v {b{v)p) = \d 2 v (a(v)p) 

for some functions &(•) and a(-), which is a second order equation. 

5. Conclusions 

The idea of interrogating legacy codes as a control system is very attractive and to some 
extend, has already been commonly used in some disciplines. For example, chemists use 
packages such as CHARMM and AMBER as legacy codes to perform optimization tasks, 
e.g. to find free energy surfaces and minimum free energy paths. Optimization techniques 
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such as the Nelder-Mead algorithm that use only function values (not the derivatives) were 
designed with this kind of problems in mind. One purpose of the work of Keller et al. is 
to extend bifurcation analysis tools to systems that are defined by legacy codes [26]. The 
"equation-free" approach attempts to extend such practices to another direction, namely 
the modeling of macroscale spatial/temporal dynamics of systems defined by microscopic 
models, in the form of legacy codes. While this seems very attractive, the set of tools pro- 
posed under this umbrella are quite far from being sufficient for reaching this objective. We 
have discussed some of the technical difficulties in this note. This discussion is certainly 
not exhaustive. It is only meant to be illustrative. 
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